#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Sun Nov  5 18:09:47 2023

@author: liqingsimac
"""
import numpy as np
from scipy.stats import zscore

c=np.loadtxt('ti12_2.txt')
n=c.shape[0]
d=zscore(c,ddof=1)          #数据标准化
r=np.corrcoef(d.T)          #求相关系数矩阵

val,vec=np.linalg.eig(r)
cs=np.cumsum(val)      #求特征值的累加和

rate=val/cs[-1]        #求贡献率
Trate=np.cumsum(rate)  #求累积贡献率
print('特征值为：',val)
print('贡献率为：',rate)
print('累积贡献率为：',np.round(Trate,4))

'''
######################################################
## 手工计算：
A0=vec*np.sqrt(val)
A=A0[:,0:3]
P=      #求出正交矩阵
A1=A@P  #因子旋转，使得 A1的列向量的分量的方差最大
A1c=np.sum(A1**2,axis=0)  #每列求和，求信息贡献
A1r=np.sum(A1**2,axis=1)  #每行求和，求共同度
print('载荷矩阵为：\n',A1)
print('对x的贡献为：',A1c)
print('共同度为：',A1r)
'''

######################################################
## 直接调用factor_analyzer包：
from factor_analyzer import FactorAnalyzer as FA
fa = FA(3,rotation='varimax')  #构建模型
fa.fit(d)         #求解方差最大的模型
A2=fa.loadings_    #提取载荷矩阵
print('载荷矩阵为：\n',A2)

gx=np.sum(A2**2, axis=0)    #计算信息贡献
s2=np.sum(A2**2, axis=1)    #计算共同度
ss=np.linalg.inv(np.diag(1-s2))
f=ss@A2@np.linalg.inv(A2.T@ss@A2)  #计算因子得分函数系数
print('共同度为：',np.round(s2,4))
print('各因子的方差贡献：',np.round(gx,4))

